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PREFACE 


The Research Achievements Reviews document re- 
search accomplished by the laboratories of Marshall 
Space Flight Center. Each review covers one or two 
fields of research and attempts to present the results 
in a form readily useable by specialists, system en- 
gineers, and program managers. 

Reviews of this fourth series are designated Volume 
IV and will span the period from May 1970 through 
May 1972. 

In accordance with NASA policy the International Sys- 
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are used in this publication. 


The papers in this report warm presented May 28, 1970 


William G. Johnson 
Director 
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STATISTICAL DETERMINATION OF SAFETY FACTOR REQUIREMENTS 
FOR UNTESTED STRUCTURES FROM SATURN V TEST DATA 

By 

Jerrell M. Thomas 


INTRODUCTION 

A traditional method of assuring structural relia- 
bility in aircraft, missiles, spacecraft, and launch 
vehicles is the application of structural safety factors. 
These safety factors must always be carefully defined 
to meet the needs of the particular system; but for 
our purposes, it is sufficient to define the safety fac- 
tor as the ratio of the load required to produce failure 
in a structure to the maximum expected load in ser- 
vice (limit load) . Typical safety factors in the aero- 
space industry range from 1.2 to 1.5 for major 
structures. It is standard practice to test a specimen 
of every major structure to assure that the required 
safety factors do exist. 

The increasing size and complexity of aerospace 
structures, coupled with the more complex combina- 
tions of environments, cause the cost of testing and 
the influence of the testing on overall development 
schedules to continually increase. This provides 
sufficient motivation to search for ways to minimize 
the test program. Is it possible to decrease the scope 
of required testing to gain cost and schedule advan- 
tages without sacrificing the all important system 
reliability? In the large systems of the immediate 
future, such as the Space Shuttle and Space Station, 
certain structures will almost certainly be sized by 
stiffness requirements from dynamic considerations 
or by long life requirements. Thus, these structures 
may have an inherently large safety factor on the 
maximum expected load. In other systems, it may be 
feasible to increase the safety factor to some desired 
level at the expense of adding weight. The purpose of 
this paper is to illustrate a possible method for deter- 
mining what safety factor is required to provide 
reliable structures without testing. 

DESCRIPTION OF METHOD 


There will' always be inaccuracies in predicting 
the failure load of any given structure. It is reason- 
able to expect that the failure load will be overpredicted 
in some cases and underpredicted in others, with most 
structures failing near the predicted value and with 


a decrease in the probability of failure at load levels 
either above or below the predicted level. As will be 
shown later, the distribution of failure frequency 
about a mean value (near the predicted value) is 
approximately normal for Saturn V structures. This 
is illustrated schematically in Figure la. Design 
changes are usually made in those structures that 
fail at less than the predicted level while those that 
fail at or above the predicted level are left unchanged. 
As a result, the failure frequency of structures that 
are actually used in flight would appear somewhat as 
illustrated in Figure lb. This odd distribution would 

LLI 

as 



PERCENT OF PREDICTED FAILURE LOAD 
a. Structures as tested 



PERCENT OF PREDICTED FAILURE LOAD 
b. Structures after testing and repair 

Figure 1. Schematic of failure frequency 
distributions. 
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then make a statistical prediction of probability of 
failure in flight very difficult, even if other 
important variables such as manufacturing dif- 
ferences and expected load uncertainties are 
ignored. 

The concern of the present analysis is to accom- 
plish the somewhat easier task of determining the 
safety factor for untested structures that will give 
essentially the same probability of success as that 
obtained with tested structures. This is accomplished 
by first assuming that the Saturn V structure is 
representative of the structural article under con- 
sideration. Statistical properties of the Saturn V test 
can be determined, and a normal distribution curve 
similar to Figure la can be constructed. Once this 
is done, a percentage of the predicted failure load 
can be determined for which a selected probability 
exists that the failure load of the article under con- 
sideration will exceed. From this percentage, the 
required safety factor can be determined to assure 
that the structural article will not fail at less than 
the ultimate load. 


TABLE 1. SATURN V TEST DATA 


Data 

Point 

Number 

Test 

Safety 

Factor 

Percent 

of 

Predicted 

Tested 

to 

Failure 

1 

1.40 

100 

No 

2 

1.58 

113 

Yes 

3 

1. 12 

80 

Yes 

4 

1 

1 

i 

1.65 

1 

1 

118 

1 

1 

1 

Yes 

1 

1 

1 

1 

1 

1 

50 

l 

1 

i 

1 

1.41 

1 

1 

1 

101 

1 

1 

1 

YeS 


Most of the Saturn V structures are made from 
high strength aluminum materials. The const "uction 
is a mixture of riveted skin-stringer, integrally 
stiffened, and sandwich articles. 


SATURN V TEST DATA ANALYSIS 


Assumptions and Limitations 

It is assumed that all structural components can 
be grouped statistically and that all tests were valid 
with respect to loads and environments. If a test 
specimen was not tested to failure, the maximum 
load sustained is taken as the strength of the speci- 
men. If one test specimen was tested to a series of 
ultimate tests, the maximum load sustained in any 
test is taken as the strength. 

Data Analyzed 

A sample of the Saturn V test data that were 
analyzed is shown in Table 1. Basically, each 
data point represents one major Saturn V struc- 
tural assembly such as a tank, interstage, or 
thrust structure. In some cases structural arti- 
cles were subjected to more than one load 
condition. If these load conditions were very sim- 
ilar to e**jh other, the maximum load sustained in 
any test was taken as the specimen strength. 

In other cases a single structural article was sub- 
jected ti more than one test, and the test 
conditions were entirely different. In these cases 
more than one data point was gained from the single 
structural assembly. 


The second column in Table 1 shows the safety 
factor attained by each article in the static structural 
test. This factor represents the ratio of the load 
attained in test to the limit load. The third column, 
which is the actual data analyzed, represents the 
percentage of the predicted capability load attained 
by the article. For purposes of this analysis, the 
predicted capability was taken as the limit load times 
the safety factor requirement for the article. An 
indication as to whether or not the article was tested 
to failure is given in the last column. Slightly more 
than half of the structures tested were tested to 
failure. The fact that some structures were not tested 
to failure produces an undesirable bias in the data, 
but this is not believed to be severe because many of 
the articles that did not fail were very near failure 
based on strain gage instrumentation. 

Statistical Representation of Data 

A statistical representation of the data is shown 
in Figure 2. The normal distribution was calculated 
tty standard statistical methods. The histogram is 
shown to give a visual comparison of the actual data 
to the normal distribution curve. The data have a 
mean value of 0. 99 and a standard deviation of 0. 142. 
A more detailed representation of the statistical 
properties is shown in Figure 3. An example of. 
conversion of the abscissa to any desired safety factor 
is shown on the (B) scale in Figure 3. 





PERCENT OF PREDICTED FAILURE LOAD (A) 

PERCENT OF EXPECTED FLIGHT LOAD (EXAMPLE. SAFETY FACTOR * 1.4) (B) 


Figure 3. Predicted failure distribution of untefte taturn V type structures. 
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Analysis of Data 

Let p represent the mean of the test data. Then, 
the percentage i- the expected failure load, which is 
Z standard deviations less than the mean value, is 
p - Zct, where a is the standard deviation. Thus, for 
untested structures, the ratio of the expected failure 
load to the actual failure load for a Z a probability is 
given by 

r = 100 

p - Zo 

Then any desired value of failure load must be multi- 
plied by this ratio to obtain a design load that gives 
a Za probability of no failure at the desired failure 
load. If the safety factor of tested structures is 
represented by FSj, the equL alent safety factor, 

FS , for untested structures would be given by 

FS = R x FS. 
u t 

In Figure 4, the required safety factor, FS^, 

for untested structures has been plotted versus die 
probability of no failure, Za, for various values 
of FS ( as a parameter. It should be noted that the 

FS t - 1. 00 .ilot (represented by the 100 percent 

limit load curve) gives the probability of failure 
for untested Saturn V type structures at limit load. 
Figure 5 is a crossplot of Figure 4 from which the 
equivalence between safety factors for tested and 
untested structures can be determined directly for 
cither a. 2a or Za probability of so failure. 


CONCLUSIONS 


Although the data presented should not be inter- 
preted as absolute because of die assumptions made 
and the limited test data available, they seem to 
yield a rational guide for choosing safety factors for 
untested structures. For structures sized by stiff- 
ness considerations, life requirements, or minimum 
gage restraints or for structures for which weight 
is not of utmost importance, a no-test philosophy with 
larger analytical safety factors to provide equivalent 
reliability is worthy of careful consideration. 



U 1J Zf M U U 

OESMN SAFETY FACTOR FOR UNTESTED STRUCTURE 


Figure 4 Failure probability prediction. 



SAFETY FACTOR OF TESTED STRUCTURE 


Figure 5. Equivalent safety factors. 
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MODAL ANALYSIS OF SIlUCTURIS IY AN ITIRATIYI 
RAYLEIGH -RITZ TICNNIQUI 

By 

John R. Admire 


SUMMARY 


A method of modal analysis of structures is 
presented. The method is based on repeated applica- 
tions of the well known Rayleigh-Ritz technique using 
improved coordinate vectors for each iteration. This 
procedure is convergent, since each application of the 
Rayleigh-Ritz technique yields a better approximation 
of die modal properties of die structure. Numerical 
examples are used to demonstrate the rapid conver- 
gence of the method. 

INTRODUCTION 


Dynamic response analyses are important in the 
design of elastic structures subjected to time- 
dependent applied forces, because the loads associated 
with these responses can exceed the loads from all 
other sources. The dynamic response of a structure 
is determined by solving its equations of motion. The 
present practice in die aerospace industry is to simu- 
late the actual structure with a discrete coordinate 
model to obtain the equations of motion. Tim discrete 
coordinate equations of motion that result are coupled 
second-order differential equations. Since coupled 
differential equations are rather difficult to solve, they 
are transformed to a set of uncoupled differential 
equr Tons. The transformation is accomplished using 
the mode shapes from a modal analysis, and the 
uncoupled equations are called modal equations. Thus, 
the most important use of iata from a modal analysis 
is to simplify the dynamic response analysis. 

The modal analysis of a discrete coordinate model 
of a structure having n degrees of freedom requires 
the solution of an nth-order matrix eigenvalue problem. 
A digital computer is usually used to solve the eigen- 
value problem. The three most important items to 
consider when choosing a method to solve the eigen- 
value problem on a computer are accuracy of the 
solution, computation time '■equired, and computer 
storage required. There are many algorithms 
available for direct solution of the matrix eigenvalue 
problem [1]; however, as the number of degrees of 


freedo. » becomes large, the computation time and 
computer storage becomes excessive for most of 
these methods. 

The technique most often used in the aerospace 
industry to reduce the computation time and computer 
storage to acceptable levels for modal analysis is 
the Rayleigh-Ritz method [2] . This method permits 
an approximation of the first m modes of the structure 
to be obtained by solving an mth-order eigenvalue 
problem rather than the original nth-order eigenvalue 
problem. Thus, when the number of modes desired 
is significantly less than the number of degrees of 
freedom of the model, a great deal of computer time 
and computer storage can be saved by using this 
method. The reduction in the order of the eigenvalue 
problem Is accomplished by approximating the dis- 
placements of the discrete coordinates by a linear 
sinn of m Rayleigh-Ritz coordinate vectors. 

The theory oo which the Rayleigh-Ritz method is 
based places only minor restrictions on the coordi- 
nate vectors that can be used. Thus, an analyst has 
great freedom in the choice of the coordinate vectors. 
When the coordinate vectors are derived from the 
modal properties of the components of the structure, 
the technique is called "modal coupling" or "compo- 
nent mode synthesis". The most well known method 
based on component modal properties is due to 
Hurty [3], 

The inherent problem in using the Rayleigh-Ritz 
method is the accuracy of the approximate solutions. 
That is, in general it is necessary to use more 
Rayleigh -Ritz coordinate vectors than die number of 
modes desired in order to obtain sufficient accuracy. 
This procedure tends to defeat the purpose of the 
Rayleigh-Ritz method since the order of the eigenvalue 
problem is equal to die number of coordinate vectors 
rather than the number of modes desired. 

The accuracy problem associated with the 
Rayleigh-Ritz method can be eliminated by applying 
the method repeatedly and using improved coordinate 
vectors each time. This procedure will converge to the 
exact modes since each application of the Rayleigh-Ritz 
method will yield a closer approximation of the exact 
modes of die structure. 


5 



JOHN R ADMIRE 


Bajan, Feng, and Jaszlics |4] have suggested 
an iterative method and applied it using coordinate 
vectors derived from component modal properties. 
The iterative method presented here appears to con- 
verge more rapidly than the method of Bajan, Feng, 
and Jaszlics, at least for the example in their paper. 


ANALYTICAL „ rVELOPMENT 

Consider a discrete coordinate model of a con- 
strained structure having n degrees of freedom. The 
free, undamped equations of motion of such a model 
expressed in matrix notation are: 

+ IK]{X} = {0} , (1) 

where {JT} is a vector of discrete coordinate displace- 
ments, (M] is the mass matrix, andlK] is the stiffness 
matrix. The time dependence in equation (1) can be 

removed by the transformation {X} =e 1Ut {x}, since 
its solutions are harmonic in time. Thus, equation (1) 
becomes 

(iKl - [M]){X} ={0} . (2) 

Equation (2) is recognized as a matrix eigenvalue prob- 
lem of order n whose eigenvectors, l YJ , are the mode 
shapes and whose eigenvalues, ['ail. ] , are the squares 
of the frequencies. 

The first step in using the Rayleigh-Ritz technique 
[2] to obtain approximate solutions of equation (2) con- 
sists of selecting a complete sequence of coordinate 
vectors 

{Z}„ {Z},,. . .{Z} m ,...{Z} n , (3) 

which are linearly independent and satisfy the kinematic 
boundary conditions. The displacement vector {X} is 
then expressed as a linear sum of the first m coordinate 
vectors 

(X}=Z<J U {Z}. . (4) 

j=‘ 

In matrix notation, equation (4) can be expressed as 

{X} = [Zl{q} . (5) 

It is known from energy considerations [2] that the 
function 

u* = {X} T IK]{X}/{X} T 1MJ{X> (6) 

must be stationary. The equations that will yield w 1 


stationary are obtained by expressing equation (6) in 
termsof{q}, differentiating it with respect to the 
Rayleigh-Ritz coordinates, q , and setting the 
result equal to 0; i. e. , 1 

a rjjf 

“ = 0 , j = 1 to m. (7) 

#, J 

Performing the operations indicated by equation (7) 
yields 

([Z] X lKJlZj - w : lZJ T lMJ|Z]){q} = (0>. (8) 

Equation (8) is a matrix eigenvalue problem of order 
m whose eigenvectors are [y* ] and whose bigen- 
values are . The approximate mode shapes 
are obtained from the relationship 

IZJ = JZ J|Y* J . (9) 

The accuracy of the mode shapes and frequencies 
obtained from die Rayleigh-Ritz method depends entirely 
on the coordinate vec'ors, [ Z j. In general, there are 
no procedures available for choosing a set of coordinate 
vectors that will yield exact solutions. However, it is 
possible to stale characteristics such a set must pos- 
sess. Since the mode shapes derived from the method 
are a linear combination of the coordinate vectors, it is 
clear that exact modes would result if each coordinate 
vector was a linear combination of the first m mode 
shapes of the structure. Another way of stating this is 
that exact results would be obtained for the first m modes 
of the structure if the coordinate vectors did not contain 
any contribution of modes higher than m. 

Consider an iterative procedure to eliminate the 
accuracy problem associated with the Rayleigh-Ritz 
technique. If the Rayleigh-Ritz technique was repeat- 
edly applied and if each time it was applied the 
coordinate vectors used would yield more accurate 
results than the previous set of vectors, it is clear 
that the iteration procedure would converge to the 
exact solution. The only problem with the iteration 
procedure is generating an improved set of coordi- 
nate vectors; that is, a set of vectors which will 
yield more accurate results than the previous set of 
coordinate vectors. Before an improved set of coordi- 
nate vectors can be generated, it is necessary to estab- 
lish the characteristics such a set must possess. 

From the discussion of the characteristics that 
a perfect set of coordinate vectors must possess, it is 
clear that an improved set of vectors is a set in which 
each vector contains less contribution from the higher 
modes than the corresponding vector in the previous 
set. Since the approximate mode shape's obtained by 
the Rayleigh-Ritz technique are a linear sum of the 
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coordinate vectors, it is possible to state the charac- 
teristics cf an improved set of coordinate vectors in 
terms of a set of approximate mode shapes; that is, 
a set of coordinate vectors will be an improved set 
with respect to the previous set of vectors, if each 
vector contains less contribution from the higher 
modes than the corresponding mode shape derived 
from the previous set of coordinate vectors. Thus, 
an improved set of coordinate vectors can be gener- 
ated by suppressing the contribution of the higher 
modes in the approximate mode shapes. The proce- 
dure for supp res sine the contribution ot higher 
monies is well known; in fact, it is the basis of the 
matrix iteration method (5] of modal analysis. 

Consider generating a vector {f} from a vector 
{s} by the following equation: 

{f> = IK] - ' lM]{s} . (10) 

The vector { s} can be expressed as a linear sum of 
the mode shapes of the structure 

n 

{s} = z (y}=b • (11) 

i=i 1 1 

Combining equations (10) and (11) yields 
n 

{f} = £ IK)' 1 ]M){y} b (12) 

j=l 1 1 

Since {y} . and m. must satisfy equation (2) , the 
following relationship exists: 

!Kj-‘ [M]{y} = ^{y} (13) 

j J 

Therefore, equation (12) can be written as 

W'Sjw. • < 14 > 

j=l j 

Comparing equations (11) and (14), it is evident 
that the contribution of the higher modes in {f} is 
less than in {s}, since the magnitude of u>* increases 

with j. Thus, an improved set of coordinate vectors 
can be obtained from the relationship 

IZ1 = (Kj-UMjm. . (15) 

It is not advisable to use equation (IS) directly to 
compute [Z] since ( K J — 1 is, in general, a full matrix 
ami requires a great deal of storage. Thus, equation 
(15) is replaced by two equations and [Z] is deter- 
mined by solving a set of simultaneous equations; 
that is. 


[F] = (M) IY] (16) 

and 

C F J = IK) |Z] , (17) 

Equations (16) and (17) can be solved efficiently for 
structures whose mass and stiffness matrices are 
banded by using methods based on Cholesky's square- 
root technique (6). 

The equations necessary in applying the iterative 
Rayleigh-Ritz technique have been developed. The 
essential steps ot the technique can be summarized 
as follows: 

1 . Select a set ot coordinate vectors [Z J . 

2. Form the eigenvalue problem expressed by 
equation (8). 

3. Solve equation (8) to yield the eigenvectors 
(Y*j and eigenvalues I'oilj . 

4. Compute a set of approximate mode shapes 
by equation (9). 

5. Compute an improved set of coordinate 
vectors [Z] by solving equations (16) and (17). 

6. Return to step 2. 

The iteration is halted when the mode shapes and 
frequencies have attained the desired accuracy. 

NUMERICAL RESULTS 


The preceding analysis was programmed for an 
IBM 7094 digital computer. Included in the program 
was a provision for computing static deflections 
resulting from unit loads applied at selected discrete 
coordinates. These static deflections were used as 
the initial coordinate vectors for the analysis. 

Modal data were generated for the cantilever 
frame structure shown in Figure 1. The structure is 
composed of 50 identical members, each of which can 
deform axially and in bending. The axial stiffness, 
AE, bending stiffness, El, and the mass density, p, 
ot the members are also shown in Figure 1 . A dis- 
crete coordinate model of the structure was generated 
using 3 degrees of freedom at each junction of the 
members, which resulted in a total of 90 degrees of 
freedom. The degrees of freedom used at the junc- 
tions were a rotation, a horizontal displacement, and 
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a vertical displacement. Modal data for this model modes, and the fourth mode is entirely an axial 
were computed using a standard eigenvalue routine mode, 

to verify the results obtained from the iterative 

Rayleigh-Ritz technique. The first three modes of Figure 2 shows a plot of the square of the 

vibration of this structure are primarily bending frequency versusthe number of iterations. Also 
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THE DISCRETE COORDINATE MOOCL OF THIS STRUCTURE HAS 3 DEGREES OF FREEDOM AT EACH 
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Figure 2. Plot of the square of the frequency versus the number of iterations. 


Figure 1. Frame structure. 
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shown in this figure is a sketch snowing the structure 
with three arrows on it indicating the locations at 
which unit loads were applied to obtain the coordi- 
nate vectors. Since three coordinate vectors were 
used, three frequencies were obtained. They are 
shown plotted with a subscript indicating the mode 
number. The values of the frequencies shown for the 
Oth iteration are the frequencies computed from the 
initial set of coordinate vectors. All three frequen- 
cies converged to a constant value after three 
iterations with the third frequency requiring the 
most iterations since it was initially in error the 
most. From Figure 2 one could conclude that if 
the initial frequencies are not greatly in error, the 
technique produces rapid convergence. The next 
question that must be answered is: if the initial 
frequencies ~re greatly in error, will the technique 
still converge to the correct solution rapidly? 

The purpose of Figure 3 is to show the rate of 
convergence of the frequencies for the case where 


the initial frequencies are greatly in error. This 
figure is similar to Figure 2 except that a different 
set of initial coordinate vectors was used. The 
locations of the unit loads used to generate the 
coordinate vectors are shown as arrows in the 
figure. These initial coordinate vectors yielded 
frequencies that were as much as two orders of 
magnitude in error. The technique still yielded 
accurate solutions after only four iterations. Thus, 
one could conclude from Figures 2 and 3 that the 
number of iterations required to obtain accurate 
solutions does not depend greatly on the initial set 
of coordinate vectors. 

Thus far, only the convergence of the frequen- 
cies has been discussed The rate of convergence 
of the mode shapes is also important. Figure 4 
shows the rate of convergence of the mode shape 
corresponding to the third frequency shown in 
Figure 3. This mode was selected since its 
frequency converged the slowest of the three shown 
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Figure 3. Rate of convergence of frequencies where the initial frequencies are greatly in error. 
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Figure 5 shows what can happen if the coordi- 
nate vectors are generated arbitrarily. Again, the 
squares of the frequencies are plotted versus the 
iteration number, and the static loads used to 
generate the coordinate vectors are shown on the 
sketch. In the previous examples, all of the static 
loads used to generate the coordinate vectors 
caused the structure to deform at least to some 
extent in bending. That is not the case in this 
example. The horizontal force shown m Figure 5 
causes no deformation normal to the centerline of 
the structure. Thus, the coordinate vector obtained 
from the horizontal force gives no approximation of 
the first three modes since they are bending modes. 
Figure 5 shows that the frequencies converged to a 
constant value rapidly, but the value to which the 
third frequency converged is not the true third 
frequency but is the fourth frequency of the structure. 
As mentioned previously, the fourth mode is an 
axial mode. Thus, the coordinate vectors chosen 
must approximate, at least to some extent, the 
type of modes that are being computed. After eight 
iterations, roundoff errors caused this mode to 
pick up a small amount of bending. When this 
happened, the iteration technique amplified this 
bending and suppressed the axial deformation, and 
it quickly converged to the true third frequency of 
the strc ture. Thus, some care must be taken 
in the coordinate vectors. 



Figure 4. Mode shape of third mode. 


in Figure 3. Figure 4 shows a plot of the centerline 
mode shape for iterations from zero to five. The 
mode shapes for the fourth and fifth iterations are 
almost identical. This indicates that the mode con- 
verged to a very accurate value in four iterations. 
Thus, from the numerical examples shown, it is evi- 
dent that the iteration techniques will yield accurate 
mode shapes and frequencies after very few iterations. 


CONCLUSIONS 


Based on the numerical examples presented, it 
appears the iterative Rayleigh-Ritz technique is an 
efficient means of computing modal data for struc- 
tures having banded mass and stiffness matrices if 
reasonable care is used in selecting the coordinate 
vectors'. 

Initial coordinate vectors obtained from static 
loads were only used to demonstrate the iteration 
technique, and it should not be concluded that they 
are the best ones to be used. It is most likely that 
initial coordinate vectors derived from the modal 
properties of the components of the structure would 
be a better choice. 

Even though the technique as presented 
applies only to constrained structures, it can be 
applied to unconstrained structures with only minor 
modifications. 
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ITERATION NUMBER 

Figure 5. Plot of the square of the frequency versus the number of iterations with the coordinate vectors 

generated arbitrarily. 


The most logical extension of this work would 
be to apply the iterative technique using initial 
coordinate vectors derived from the modal properties 


of the structure’s components and generate the 
improved coordinate vectors using substructuring 
techniques. 
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FLOW INDUCED VIBRATIONS IN METAL BELLOWS 

By 

H. J. Bandgren* 


SUMMARY 

This paper summarizes the results of an experi- 
mental and analytical study of flow induced vibration 
of metal bellows. The results of numerous labora- 
tory experiments have confirmed that vortex shedding 
from the tips of the bellows convolutions is the fluid 
excitation mechanism. When the frequency of the 
vortex shedding coincides with one of the natural 
longitudinal resonant frequencies of the free bellows 
or braided hose, a strong bellows vibration may 
exist. 

A semiempirical method of estimating the flow 
induced dynamic stresses in a bellows for a given 
set of flow conditions has been formulated. Design 
methods for eliminating or minimizing the effects of 
flow induced vibration in bellows systems have been 
developed and experimentally verified. 

The criteria presented in this paper were used 
to assess all of the metal bellows in the Saturn/ 
Apollo launch vehicle system to determine if flow 
induced vibration problems existed. Corrective 
action was implemented where necessary by the 
installation of flow liners or the addition of damping 
to the bellows system. 


INTRODUCTION 


During the flight of the second Saturn V space 
vehicle, AS-502, one of tV 5 -2 engines on the second 
stage malfunctioned causir;' die engine to shut down 
early. During the same flight the J-2 engine on the 
third stage failed to restart in earth orbit. More 
recently, during static firings of the S-IVB flight 
stages for AS-508 and AS-509, the LH 2 propellant 
delivery lines developed leaks. All of these anoma- 
lies were determined to be the result of flow induced 
vibration fatigue failures of the metal bellows. 

This paper describes the results of an analytical 
and experimental study of flow induced vibrations in 


metal bellows. A number of significant findings have 
been made throughout this study and are discussed in 
detail in the paper; these are listed below: 

1. The fluid-elastic mechanism causing bellows 
flow excitation (vortex shedding) has been observed 
in the laboratory and is discussed. 

2. A semiempirical method for estimating the 
flow induced dynamic stresses in a bellows for a 
given set of flow conditions is described. 

3. Methods of eliminating or minimizing flow 
induced dynamic stresses are presented. Internal 
flow liners and various bellows external damping 
devices have been tried as a means of suppressing 
the bellows flow induced vibrations. 

VORTEX SHEDDING FROM BELLOWS 
CONVOLUTIONS, THE FLOW EXCITATION 
MECHANISM 

It is documented in the literature that when there 
is fluid flow over a cylinder, as shown in Figure 1, 
with a sufficiently high Reynolds number, a succession 
of alternately shed vortices will occur in the near 
wake region of the cylinder. Over a Reynolds num- 
ber range of from 300 to 300 000, the frequency of 
the vortex shedding is related to the Strouhal number 
(S), flow velocity (V), and the diameter (D) of the 
cylinder as shown in Figure 1. This same phenomenon 
has been shown to occur when fluid flows through a 
bellows with open convolutes, as depicted at the 
bottom of Figure 1. The length quantity is a, the 
width of the convolute. 

When the vortex shedding frequency of the bellows 
coincides with one of the longitudinal resonant fre- 
quencies of the bellows, a strong bellows vibration 
may result. To verify this, the response characteris- 
tics of numerous bellows were studied in the flow 
facility. Figure 2 shows the flow induced dynamic 
stress at the tip of the convolute for a typical free 
bellows. This is a single ply bellows with an inside 
diameter of 5. 08 cm (2. 0 in. ) and a convolute 
width of 0. 3175 cm (0. 125 in. ) . A graph of root 


*This study was performed in conjunction with Southwest Research Institute under contract NAS8-21133 
with Dr. C. R. Gerlach as the principal investigator. 


13 



H. J. BANDCREN 



« = SV/D = 0.2 Y/D 

VORTEX SHEDDING FROM CYLINDER 



f = SV/(T = 0.2 v/<r 

VORTEX SHEDDING FROM BELLOWS CONVOLUTIONS 


FLOW RATE - liters/sec 



Figure 1. Vortex shedding process. 


Figure 2. Bel) jws stress as a function of flow rate 
for four ilow excited modes of vibration. 


mean square stress at the tip of the end convolution 
is plotted against flow rate in liters/sec (gallons 
per minute are also shown) . As the flow rate 
increased, a successka of bellows modes were 
excited. Before the bellows was installed in the 
flow facility, it was filled with water, capped off, 
and fixtured on an electrodynamic shaker to identify 
its longitudinal natural frequencies. These were 
400 Hz for the first mode, 730 Hz for the second 
mode, 1130 Hz for the third mode, and 14F0 Hz for 
the fourth mode. The significance here is that the 
flow rate corresponding to each of the stress build- 
ups shown in Figure 2 resulted in a particular vortex 
shedding frequency that coincided with one of the 
longitudin: 1 ' resonant frequencies of the bellows 
when a Sirouhal number of approximately 0. 2 was 
used in the calculation for the vortex shedding 
frequency. Also, the dynamic stress increased 
as the mode number increased, because the vortex 
shedding forces are proportional to the free-stream 
stagnation pressure, % pV 2 (where p is the fluid 
mass density) . 


The type of longitudinal resonant modes that 
are excited in a free bellows are shown in Figure 3. 
These resonant modes are sometimes referred to 
as the accordion modes. For the first four natural 
modes of the previous bellows, axial displacement 
i« plotted on the ordinate as a function of position 
along the bellows on the abscissa. The first 
natural mode is 400 Hz. There are nodes at the 
ends of the bellows, and all of the convolutions are 
moving in phase with the maximum motion at the 
center of the bellows. The second mode is 780 Hz. 
The nodes are at the center and end points of the 
bellows, and the maximum motion is at the quarter 
points 180 degrees out of phase. See Figure 3 for 
a description ot ihe third and fourth modes. 

A bellows vibrating in one of its accordion modes 
is better represented by a lumped parameter system 
than by a continuous system. The mechanical model 
shown at the bottom of Figure 3 was found to give 
quite adenuate predictions of the longitudinal bellows 
frequencies. The model consists of a series of equal 
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Figure 3. Observed bellows longitudinal vibration 
mode shapes. 


springs (k) and masses (m) fixed at the ends. The 
individual springs can be determined from the axial 
spring rate of the bellows: the masses can be deter- 
mined from the total mass of the bellows convolutions 
plus the mass of the fluid trapped in the root of the 
convolutes. The mass and flexibility matrices can 
be determined from the model, and any standard 
eigemunction routine will provide the res >nant 
frequencies and mode shapes. 

For a braided hose the type of modes that are 
excited by the vortex shedding forces are shown in 
Figure 4. The adjacent crowns of the convolutes are 
fixed by the wire braid, and a simple cantilever mode 
of the individual roots of the convolutes is excited. 
These are sometimes referred to as reed modes, and 
the spring mass system is a simple one as shown at 
the bottom of Figure 4. 

Experimental stroboscopic observations have 
been used to verify that these are the longitudinal 
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Figure 4. Metal hose flow excited vibration mode. 


bellows modes that are excited by the vortex shedding 
forces. To experimentally verify that vertex shedding 
is the mechanism that is exciting the longitudinal 
modes of the bellows, a model consisting of a two- 
dimensional clear plastic channel containing a short 
section of convoluted metal was constructed, as 
shown in Figure 5, and placed in the flow facility. 

As in an actual bellows, the convolution tips are 
exposed to the fluid flow. When water is passed 
through the channel at t l 3 proper flow velocity, the 
convolutions are flow excited. By injecting ink 
'pstream from the convoluted segment and by slowing 
down the motion with the aid of a strobe light, the 
vortex shedding process was readily visualized, and 
motion pictures were taken. 

Figu.e 6 shows the sequence of fluid and convolu- 
tion events that were observed in a frame-by-frame 
examination of the motion pictures. Note that the 
mode of vibration of the segment is one where each 
convolution moves out of phase with the adjacent 
convolutions. The vortex shedding process on the 
vibrating convolution segment, as shown in Figure 0, 
occurs as follows: 

1. Position I — A large vortex c has formed 
between convolutions 1 and 2. A small vortex b is 
beginning to term on the downstream side of convolu- 
tion 2. A large vortex a is moving across the tip 
of convolution 3; its origin will become apparent 
later. 


15 



7 

INK INJECTION PORT 




SLIDING 


BLOCK & ROD 


2.54 tm (1.0 is.) WIDE STRIP CUT FROM 
20.32 ca {1.0 ia.) IILIOWS 


Figure 5. Two-dimensional bellows flow visualization model. 


2. Position II — Vortex r has been pushed 
downstream by the pinching action of convolutions 1 
and 2. Vortex b is gaining in strength on the down- 
stream side of convolution 2. Vortex a is propagating 
downstream. 

3. Position III — The space between convolu- 
tions 1 and 2 has opened up causing the fluid to be 
drawn in, and vortex d is beginning to form on the 
downstream side of convolution 1. Vortex c has 
been swept over the- top of convolution 2 with the 
timing being about right to combine with vortex b, 
forming the larger vortex b + c. 

4. Position IV — Vortex d is gaining in strength 
on the downstream side of convolution i. Vortex 

b + r : is being pushed out into the fluid stream by the 
"pinching" action of convolutions 2 and a. 

5. Position V — This is the same as position I 
completing one cycle of the vortex shedding process. 
Note that vortex a in position 1 is the combination of 
two vortices, as is b + c in position V. 


SEMIEMPIRICAL METHOD FOR ESTIMATING 
FLOW INDUCED DYNAMIC STRESSES 
IN METAL BELLOWS 


The previous discussion has been directed to the 
problem of when and how a bellows car be flow 
excited. The coincidence of the longitudinal resonant 
frequencies of a bellows with the vortex shedding 
frequency does not indicate the severeness of the 
vibration. During this study, an attempt has been 
made to predict vibration amplitudes and correspond- 
ing stress levels in a bello ws for a given set of flow 
conditions. 

Figure 7 shows the fluid and structural dynamics 
involved in bellows flow excitation. The process of 
periodic vortex formation and shedding causes a 
corresponding periodic pressure to be exertsc on 
the bellows convolutions. The amplitude of this 
alternating pressure is proportional to the free-stream 
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Figure 6. Sequence of coupled fluid-convolution 
events observed with two-dimensional bellows 
flow visualization model. 

stagnation pressure % pV 2 . So far as the bellows 
structure is concerned, the effect of this alternating 
pressure may be considered as a net force applied at 
the tip of each convolution. This force is called the 
vortex shedding force. The semiempirical equations 
for the vortex force, the convolution displacement, 
and the resulting stress at the tip of the convolution 
are shown in Figure 7. 

Figure 8 shows the model that was used to 

determine the value of the vortex force coefficient 

C_. A single be' lows convolution was simulated 
r 

in this test by a flexible steel ring clamped in a 
special housing between a pair of exciter coils. 

A displacement probe was built into the apparatus to 
allow ring vibration amplitudes to be monitored. 

The convolution vibration test model was placed in 
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the flow facility, and the vibratory response of the 
ring in water was determined with the excitation 
provided by the exciter coils. Next, the fre- 
quency response of the ring was determined with 
the fluid flow providing the excitation. Finally, a 
static force versus deflection calibration was made 
on the ring. The data from these tests were used 
as follows: 

1. Vibration amplitude as a function of fluid 
velocity was converted to an apparent force using 
the data from the force-deflection calibration. 

2. This apparent force was converted to a true 
vortex force by dividing it by the dynamic amplifica- 
tion factor Q, obtained from the forced vibration test 
in water. 


3. The vortex force coefficient C was then 

determined from the equation for the vortex force 
shown in Figure 7. 

Voitex coefficient data were obtained' in this 
manner for several convolution geometries. They 
were also obtained from numerous tests on real 
bellows. The results of these tests are summarized 
in Figure 9. The vortex force coefficient is 

plotted on the ordinate as a function of the ratio of 
pitch to convolute tip width on the abscissa. Notice 
the great reduction ia the force coefficient as the 
convolution pitch is opened. 

The next factor in the vortex force equation 

shown in Figure 7 is the elbow coefficient C . The 

£ 



Figure 9. Summary of bellows vortex force coefficient experimental data. 


18 




H. J. BANOGREN 


presence of an elbow upstream front a bellows causes 
excessive flow induced dynamic stresses because of 
the increased turbulence. There is little quantitative 
data available on this effect, but the available data 
suggest correcting the vortex shedding force upward 
by ;> factor of 2. 0 to account for the presence of an 
upstream elbow. 

A is the projected area of the convolution and 
P 

can be approximated by the product of the mean cir- 
cumference and the height of the convolute. The 
final quantity in the vortex force equation is the 
free stream stagnation pressure '/ 2 pV 2 . 

The convolution displacement resulting from 

the vortex force is shown in the second equation of 

Figure 7. The vibration mode factor C relates 

m 

displacement and force for a given resonant condition. 
The following assumptions were made in the deriva- 
tion of the vibration mode factor: 

1. The shape of the longitudinal bellows modes 
is intermediate between a linear and sinusoidal 
form. 


2. The maximum deflection for a given mode 
is equal to the static deflection times Q, the dynamic 
amplification factor. 

3. The maximum stress point for a given 
vibration mode occurs at the convolution with die 
maximum relative displacement. Generally, this 

is the end convolution for each half-wave mode shape. 

The first assumption implies that the mode shape 
over the first quarter wavelength is 

~ f [(?)»♦ *(¥)] ■ 

where x denotes the axial absolute displacement of 
a given point along the bellows defined by the axial 
position coordinate Y, N is the mode number, and 
i is the length of the bellows. 


The second assumption implies that the maximum 
absolute displacement x° is 


x 

o 


FQ 

4NK, 

A 


( 2 ) 


From the third assumption, the point of maxi- 
mum relative displacement for a given bellows, 

and a given mode of vibration, occurs at a point 
where Y - </2N„. Therefore, from equation (1) 



Combining equations (2) and (3) yields 



and by comparison with the equation for the convolu- 
tion displacement in Figure 7, the vibration mode 

factor C is 
m 



Returning now to the equation for the convolution 
displacement shown in Figure 7, the next factor is 
the vortex shedding force F that has been discussed 
and then the dynamic amplification factor Q, a 
measure erf damping. Physically, Q represents the 
ratio of the vibration amplitude at resonance to the 
vibration amplitude in the limits of zero frequency. 

It also can be expressed as the resonant frequency 
divided by the frequency bandwidth at the half-power 
points for a given resonant buildup. Because Q is so 
important in the prediction erf bellows flow induced 
dynamic stresses, a large number of tests were 
performed from which Q values could be obtained. 
Each bellows tested was installed in a special fixture 
that rigidly connected the end flanges. This was to 
insure that the vibration modes excited were the 
same as those excited by the vortex shedding forces. 
The bellows were then excited with an electrodynamic 
shaker with die output strain at the tip of the convo- 
lutes measured to obtain Q values. Some observa- 
tions noted during this testing were as follows: 

1. As the strain was increased, the damping 
increased and lowered the Q values. 

2. Increasing the number of plies in a bellows 
r reases the damping significantly. 

3. Going from a low pressure gas to a high 
pressure gas or from a light liquid to a dense liquid 
increases the damping. 


where F is the vortex shedding force, Q is the dyna- 
mic amplification factor, and K is the axial spring 
rate of the bellows. 


Some typical Q values as a function of number 
of plies and flow media are: 
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No. Plies 

Medium 

Q Range 

1 

Air or low pressure gases 

150 to 200 

1 

High pressure gases or 
lifd>t liquids 

100 to 150 

1 

Heavy liquids 

50 to 100 

2 

Gases 

40 to 50 

2 

Liquids 

30 to 40 

3 

All 

20 to 40 


The last factor in the convolution displacement 
equation is the axial spring rate of the bellows. 

The most accurate method for determiningK. is 

A 

from load deflection tests in the laboratory. Probably 
the most accurate analytical method for determining 
bellows spring rate was developed by the Battelle 
Memorial Institute. This is a computerized method 
that provides for both linear and nonlinear elastic 
deformations. Most noncomputerized techniques 
use equations s imila r to those developed by F 
Salzmann. 

The last equation in Figure 7 is for the resultant 
stress at the tip of the bellows convolution. The stress 
is related to Young's modulus (E) and ply thickness 
(t) of the bellows material, and the height (h) and 
deflection of the convolute by a geometric stress 
factor C g - Here again, the Battelle Memorial 

Institute has obtained very good agreement between 
theoretical and experimental methods using a com- 
puter solution for the stress deflection r el ations h ip. 
However, the Salzmann equations for determining 
dte stress in a bellows for a given convolution dis- 
placement compare favorably with the experimental 
test results. 

METHODS OF ELIMINATING OR MINIMIZING 
FLOW INDUCED DYNAMIC STRESSES 

The mechanism by which metal bellows are 
excited by flow induced vibration and a semiempirical 
method of estimating die resulting stresses have 


been discussed. The question now is: how can 
these dynamic stresses be eliminated or minimized? 
The installation of an internal flow liner in a bellows 
will eliminate flow induced dynamic stresses , when 
the liner design covers all of the active convolutions. 
The cone shape of most liners allows angular 
movement of the bellows but produces a constriction, 
resulting in unwanted pressure drops. 

A method for reducing the flow stresses ii a 
bellows is tc add damping to the bellows structure. 
This can be accomplished by the following methods: 

1. Increasing the number of plies of the bellows. 

2. Filling the convolutions with some visco- 
elastic material such as rubber or 3M strip-caulk. 

3. Placing a soft metal spring in each of the 
bellows convolutions. 

4. Tightly wrapping the exterior of the bellows 
with wire screen. 

Experimental results show all of these methods 
afford a considerable amount of damping. 

CONCLUSIONS 

An examination of test results of flow induced 
vibrations of numerous bellows confirmed that 
vortex shedding from the bellows convolutions is the 
excitation mechanism. When the vortex shedding 
frequency coincides with one of the natural resonant 
modes of the bellows, a strong vibration may exist. 

Tests results compare favorably with the 
semiempirical method presented in this paper for 
predicting Hie flow induced dynamic stresses. 

They also confirm that the external damping devices 
suppress flow induced vibrations of bellows. 

The results of this study were used to assess 
all of the bellows on the Saturn/Apollo launch 
vehicle system to determine flow induced vibration 
susceptibility and to provide guidelines for all 
necessary corrective action. 
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NASTRAN, NASA’S GENERAL PURPOSE STRUCTURAL 
ANALYSIS PROGRAM 

By 

R. L. McComas 


NASTRAN is a general purpose digital computer 
program for the analysis of large complex structures. 
The acronym NASTRAN is formed from NASA struc- 
tural analysis. During the annual review of NASA's 
research program in the area of structural dynamics 
in January 1964, it became apparent that there was 
considerable effort by many of the Centers to 
develop computer programs for structural analysis, 
designed to meet the particular needs of each of the 
Centers. It was suggested that perhaps a single 
program could meet all their needs. The Office of 
Advanced Research and Technology appointed a 
committee with representation from eight NASA 
Centers to study this possiblity. Tim:;, the Ad Hoc 
Group on Computer Methods in Structural Analysis 
was formed. 

After 6 months of investigation, the Ad Hoc 
Group reported to NASA Headquarters that there 
was no digital program in existence that had broad, 
uniform capabilities in the three interdependent 
disciplines of analytical mechanics, numerical 
methods, and computer programming. The Group 
did observe that there was considerable capability 
dispersed throughout the aerospace industry that 
had not been coUected into a single program. They 
also found that there was a tendency toward pro- 
prietary secrecy that inhibited any exchange of 
information. Communication was further hindered 
by the lack of compatibility between the structural 
analysis programs of any two companies. The 
Ad Hoc Group recommended that NASA sponsor an 
entirely new program aimed at bringing together 
all the best structural analysis computer techniques 
in the state-of-the-art. 

NASA Headquarters endorsed the recommenda- 
tions of the Ad Hoc Group and commissioned them 
to prepare a set of specifications. Fortunately, the 
papers from the first Wright Field Conference, 

"On Matrix Methods in Structural Analysis", were 
available to consult. The foUowing objectives of 
the specifications have now been established: 

1. Organize the program to be general-purpose. 

2. Embody a large three-dimensional struc- 
tural capability. 


3. Establish computer independence. 

4. Provide for modification without cascading 
effects. 

5. Build-in the maximum of user convenience. 

6. Document all aspects to gain maximum 
visibility. 

The contract to implement the NASTRAN 
specifications was awarded to Computer Sciences 
Corporation (CSC) with MacNeal Schwendler, 

Martin Baltimore, and later BeU Aerosystems 
Company as subcontractors. The quality of the 
NASTRAN program and its documentation is testi- 
mony to the purposefulness with which the members 
of the implementation team applied themselves. 

This team has often exceeded the state-of-the-art 
guidelines that were established. A few examples 
are the segment file allocator, the general input/ 
output module, matrix decomposition with active 
columns, the inclusion of scalar nonlinearities in 
control of system dynamics, the generality in the 
plot module, and the development of the self- 
contained "inverse power with shifts" module for 
eigenvalue extraction. The overall design of the 
program has set a new standard for general purpose 
programs of any discipline. The framework used in 
NASTRAN can be disassociated from elastic structures 
and can be applied to other disciplines, because there 
are no semantic implications in the executive opera- 
tions. The program abounds in service code, which 
threads through every step of the problem physics 
providing convenience to the analyst. 

Many policy decisions had to be made as to the 
content of NASTRAN. The total framework of the 
program was considered to be the most important; 
the executive system had to be capable of managing 
problems unbounded by core, be compact in its 
space requirements, be able to restart problems, 
be able to operate efficiently over different com- 
puters, and still be maintainable. 

It was decided that only basic finite elements 
would be included to deal with one- and two-dimensional 
elastic relationships such as beams, plates, and 
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axisymmetric shells. Economy of running time 
for the state-of-the-art elements was the 
determining factor in their selection. More 
sophisticated and new elements have appeared 
in the state-of-the-art since the inception of this 
program; these can easily be incorporated. 

The decision was made to write a single program 
in Fortran IV with some exceptional areas to be 
written in assembly language. Fortran .TV, version 
13, in its various forms on different computers, 
seems to have stabilized as a language that will be 
current for a reasonably long time. More than 
99 percent of the program is written in Fortran. 

The program is modular, so that updating is a 
matter of revamping material within a module with- 
out modifying its external appearances. 

Hopefully, the framework that has been built has 
sufficient capability to serve a sizable portion of the 
large problems that currently exist in the structural 
analysis community. It is intended that new capa- 
bility be added or outdated capability be replaced 
by augmenting or replacing modules. An entire 
chapter in the Programmer's Manual has been 
devoted to the topic of Modifications and Additions 
in anticipation of serving this particular activity. 

The traffic in the development of new and increasingly 
versatile elastic modules is expected to be the most 
active. As analysts increase their use of the program, 
their detailed needs will become better defined and 
the result will be that the traffic in "convenience . 
code" will also be expected to increase. It is desired 
that all such new features be called to the attention 
of the NASTRAN Project so that these ideas and 
routines can be disseminated to a broad audience. 

A computer program having all of the capabilities 
of NASTRAN will obviously be large, as NASTRAN 
is. In fact it is larger than the entire Univac 1108 
Exec 8 System Program at the Marshall Space Flight 
Center Computation Laboratory. 

The overall effectiveness of a general purpose 
program depends in large measure on how well thr 
available programming techniques have been employed 
in the design of its organizational and control fea- 
tures. NASTRAN has a modular separation of func- 
tional capabilities organized under an efficient 
problem-independent executive system. This 
approach is absolutely essential for any complex 
multioperation, multifile application program. 

Despite the modular separation, it is clear that 
modules cannot be completely independent since they 


are all directed toward solution of the same general 
problem. In particular they must intercommunicate 
data between themselves. The executive routine 
separates system functions from problem solution 
functions. It establishes and controls the sequence 
of module execution according to options specified 
by the user; it establishes and communicates values 
of parameters for each module; it allocates files for 
all data blocks generated during program execution 
and performs input/ output to auxiliary files for each 
module; and it maintains a full restart capability for 
restoring a program execution after either a scheduled 
or unscheduled interruption. The executive system 
is open-ended in the sense that it can accommodate 
an essentially unlimited number of functional modules, 
files and parameters. 

The functional modules am those module s that 
perform the actual calculations. They belong to one 
of four catagories; i. e. , structural modules, matrix 
operations, utility modules, and user modules. 

The structural modules are the main subprograms 
of NASTRAN. Some examples of structural modules, 
taken from the dynamic analysis, are real eigenvalue 
analysis, modal dynamic matrix assembler, transient 
dynamic response, and dynamic data recovery. The 
NASTRAN development team has made every effort 
to incorporate the very latest and best structural 
techniques to perserve accuracy and speed of solution. 

Most difficulties in numerical analysis arise in 
connection with three basic implicit operations; 
matrix decomposition (or inversion) , eigenvalue 
extraction, and integration of differential equations. 
The major difficulties that occur in the application 
of these operations to large problems are excessive 
computing time, error accumulation, and instability. 
Many methods that work well with small or moderate 
size problems are not acceptable for large problems. 
The method employed for matrix decomposition 
(triangle decomposition with active column) is 
especially important because of its extensive use as 
a base for the other two implicit operations. The 
method that is employed in the program takes maxi- 
mum advantage of matrix sparsity and bandedness. 

The latter aspect is particularly important because 
of the enormous gain in efficiency that occurs 
when banding techniques are properly employed by 
the user in setting up problems for the displacement 
method. 

In general, the solution time for a large struc- 
tural analysis of any type can be greatly reduced by 
taking full advantage of the sparsity and bandwidth 
of the matrices that describe the structural problem. 
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Other means have been used to improve efficiency fc" 
large problems. These include storing sparse 
matrices in packed form. 

The input/ output options have been designed to 
be flexible in order to accommodate the analyst. The 
plot package is very versatile allowing a choice of 
orthographic, perspective, or stereoscopic X-Y 
curve plot. All or portions of the structure can be 
plotted with varying view, line weight, and color. 

The explanation of any complex activity must be 
subdivided into phases or steps to be intelligible. In 
the case of a computer program for structural 
analysis, it is convenient to divide the total effort into 
a problem formulation phase and a problem solution 
phase. The beginning of the problem formulation 
phase occurs in the mind of the analyst. He contem- 
plates, decides what he needs to know, and constructs 
a mathematical problem whose solution, he hopes, 
will provide relevant answers to his question. The 
range of choice in mathematical problem formulation 
provided by NASTRAN is, however rich in detail, 
limited to one basic approach; namely, the use of 
finite element, structural models. The idealized 
structural model in NASTRAN consists of grid points 
to which loads are applied and at which degrees of 
freedom are defined. The grid points are connected 
by elements in the NASTRAN element library. The 
structural element is a convenient localizing concept 
for specifying many of the properties of the structure 
including material properties, mass distribution, 
and some types of applied loads. 

A wide range of elastic elements including several 
forms of beams, triangular and quadrilateral plates, 
scalar elements, and special shell and axisymmetric 
solid elements have been included in the element 
library. 

Anisotropic material properties may be employed 
in all plate elements. With modeling experience, 
there is no aerospace structure that cannot be 
accurately modeled by the NASTRAN elements. Some 
examples that are currently being solved by NASTRAN 
are the dished head assembly (Fig. 1) and the ATM 
solar panels (Fig. 2) . The idealized model of the 
dished head assembly (Fig. 3) is formed by beams 
and triangular plates. The grid mesh has been 
increased in the area of interest, the intersection of 
the two surfaces. Since the assembly is symmetric, 
only half of the structure needs to be modeled. 



Figure 1. Dished head assembly drawing. 



Figure 2. Finite element model of ATM 
solar array. 



Figure 3. Finite element model of dished head 
assembly. 

formats, as they are called in NASTRAN, for statici 
and seven for dynamic analysis. 

In tills problem , the stress and deflection? 
were desired. However, it is possible to save die 
problem solution and to restart and perform a 
different analysis at a later date, possibly a vibration 
analysis, without dupli.- ting the original effort. 


Now, all that is left to be done is to fill out the 
NASTRAN data cards (a simple although sometimes 
tedious task) and choose the analysis type, of which 
there are twelve available. There are five rigid 


The documentation for the NASTRAN Computer 
Program consists of three manuals; the Theoretical 
Manual, the User’s Manual, and the Programmer's 
Manual. 
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The Theoretical Manual is a commentary on the 
program. It is, first of all, intended to be an 
introduction to NASTRAN for all Interested persons, 
including those who will use the program and those 
whose interests are less direct. For this purpose, 
the structure and the problem-solving capabilities 
of the program are described in a narrative style. 

The manual's most impoVtant function, however, 
is to justify the analytical and numerical procedures 
that underlie the program by presenting explanations 
and derivations of sufficient depth to bt convincing. 

The selection of material has not been an easy 
task because not everyone has the same concept of 
what the word "theory" means when it is applied to 
a computer program. A broad view concerning the 
selection of materials has been adopted; i. e. , the 
formulation that will be solved, the development of 
the procedures or algorithms that will be used in the 
solution, the organization of the program, and the 
flow of data through the computer. 

The User's Manual contains all of the informa- 
tion needed to solve problems with NASTRAN. It 
includes instruction in structural modeling techniques, 
instruction in input preparation, and information to 
assist the interpretation of output. It contains 
descriptions of all input data cards, restart proce- 
dures, and diagnostic messages. 

The Programmer's Manual contains the informa- 
tion that is required for maintenance and modification 
of the program. It contains a complete description 
of the program code including the mathematical 
equations that are implemented in the Functional 
Models and describes the Executive System and the 
coding practices that have been employed. 


The NASTivAN program is currently operational 
on the IBM 7094/7040, the IBM :)G0 models 50 
through 95, the CDC 6400, 6500, and 6600, and 
the Univac 1108 Exec JI computer systems. The 
piotters that car be utilized with these systems 
are the Stromberg Carlson 4020, the Bensor. Lehner 
LTE or STE, the Electronics Associates EIA 5500, 
the Display Data DD80, and the Calcomp 765. 

Now that the initial development phase has been 
completed, there will be continuing maintenance by 
a NASTRAN Management Center at a yet undetermined 
NASA Center. The Office of Advanced Research and 
Technology will fund the maintenance and development 
tasks. In other words, the program will not become 
stagnant or die. In fact, the first year's improve- 
ments have been determined and priorities have been 
assigned. A few examples are; 

1. Substructure Partitioning 

2. Finite Element Thermal Analyzer 

3. Dummy Element 

4. Single Precision Option 

5. Training Films 

NASTRAN is the perfect tool to be used in the 
development oi the shuttle vehicle. It will, first 
of all, do the job and do it efficiently; it will allow 
the monitoring N/'SA Centers, the prime contractors, 
the subcontractors, and other interested parties to 
exchange information freely and to converse in a 
common language; and possibly of equal importance, 
with the correct management and guidelines, it 
can save untold development time and costs. 


26 



SKYLAB DOCKING RESPONSE ANALYSIS 

By 

W. B. Holland 


ABSTRACT 

The salient features of the mathematical 
formulation and digital algorithm developed by 
Martin-Marietta Corporation. Denver Division for 
the Skylab docking response analysis are presented. 
The mathematical model considers the docking .' lies 
to be in proximity such that collision is imminent. 

The coupling of the docking mechanism and the 
elastic responses of the vehicles are included. The 
Apollo probe/drogue docking system is considered. 

INTRODUCTION 


The docking maneuver and its associated 
structural responses is a primary consideration in 
the design of large flexible space structures such as 
Skylab. The Martin-Marietta Corporation, Denver 
Division under contract to Marshall Space Flight 
Center has developed the methodology and digital 


algorithm to determine the response of two elastic 
vehicles during the docking maneuver. The objective 
of this paper is to present the salient features of 
their analysis. More rigorous explanations and 
developments are contained in References 1 through 3. 

Many authors ha. 3 analyzed the docking maneu- 
ver in recent years In general, their formulation'' 
were restricted to consideration of rigid vehicles 
Their docking analysis was performed in two steps; 

(»} calculation of the force-time histories at the 
docking mechanism between the vehicles and (2) 
application of the force-time histories from step l 
to elastic response analyses of the respective bodies. 
This procedure is referred to as the "two-step 
method. " It was shown in Reference 4 that the two- 
step method can yield very conservative elastic 
responses for t highly flexible structure such a 
Skylab (Fig. l ) . More realistic responses and loads 
are attainable by applying a "one-step method", 
wherein the effects of docking elastic rather than 
rigid vehicles are included, a he Skylab Docking 
Response Analysis is a one-step method. 



Figure 1. Skylab docking configuration. 
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PHYSICAL PROBLEM 

The docking maneuver can be described as a 
controlled collision of two vehicles. The vehicles 
are considered to be initially in proximity such that 
collision is imminent. Each vehicle has a docking 
mechanism hat serves to ( 1 > dissipate the relative 
kinetic ei :gy. (2) align the vehicles, and (3) bring 
the vehicles together and thus achieve a structural 
t'3. 

The Apollo probe/drogue docking system (Fig. 2/ 
is utilized in the Skylab program. The clocking prebe 
is mounted on the command service module (CSM) , 
and a drogue is mounted on the Skylab cluster at a 
radial docking port and also at the axial docking port. 
The CSM will be referred to as the chase vehicle, 
and the cluster will be referred to as the target 
vehicle. 

The decking probe (Fig. 3) consists of a central 
body, probe head, pitch arms, tension links, hydrau- 
lic attenuators, capture latches, and support struc- 
ture. The central body contains an im-er cylinder 
that allows relative axial motion (10-inch stroke) 
with respect to an outer cylinder. The probe head 
is gknbal mounted on the inner cylinder and houses 
capture latches that engage the drogue apex. 

The pitch arms and tension links act to align 
the vehicles and prevent any jack-knifing motion. 


The attenuators are fluid-displacement type units 
and are the primary sources of energy dissipation. 
The geometry of the probe and drogue is such that 
the potential points of contact are the probe head on the 
drogue surface and the pitch arms on the drogue lip. 

ANALYTICAL DEVELOPMENT 


Curing the course of the docking simulation, 
there are constraints that result from the geometry 
of the docking mechanism. As an example, consider 
the requirement that the motion of the probe head 
be constrained to be coincident with the motion of 
the drogue apex when the probe heed is captured. 

Such conditions of constraint do not exist a priori, 
but are maintained by strong forces. It is a distinct 
advantage of the Lagrangian method 1 5 J that allows 
one to formulate the equations of motion for a system 
with constraints without knowledge of the constraint 
forces. The constraint forces are implied by 
constraint equations. Now, the form of Lagrange's 
equations with auxiliary conditions of constraint 
are cited. 

EQUATIONS CF MOTION 


Consider the dynamical system to be character- 
ized by discrete coordinates 



Figure 2. Docking system major assemblies. 
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Figure 3. Probe docking system. 


<Ji. • .Qjj ■ (!) forces required to maintain the constraint. The m 

equations of constraint coupled with the n 

The motion of the coordinates is restricted by differential equations of motion are sufficient to 

kinematical conditions. These auxiliary conditions solve for the n unknowns, q^; and for the m 

of constraint may be expressed in general form as unknowns, X . 


<t> s = <P a <«u. ii. <ii, <b, q^, t, dt) = o . (2) 

s = 1, 2, . . , m 


The Lagrange equations of motion can be 
expressed as 

d / 3T \ <1T v r 

T7 1 7 I- 7 = Q. + >. A ’ <b . 


In the above equation, T is the kinetic energy of the 
whole system. The generalized forces, Q^, are 

associated with all impressed forces both conserva- 
tive and nonconservative. The Lagrangian multiplier, 
Xg, is a measure of the microscopic violation of the 

constrant equation <(> g = 0, and is equivalent to the 


In general, it is more convenient to apply 
equations (2} and (3) in matrix form. To achieve 
a matrix form, consider a transformation 
matrix, I/Jj, where 

{«} = - (4) 

The vector {u} contains nonholonom ic velocities 
(velocity components projected on a totaling coordi- 
nate system) . 

A matrix (Pi is defined such that the kth row of 

[P] is 


(kth row of [P]) = {q} T T , (5) 

. kj 
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The auxiliary conditions of constraint (equation 2) a set of first order differential equations of the 

are expressed- in the form general form 


*s B S iA + *st = 0 - < 6 > 

k-1 

The term 0 has been added to account for rheoncmic 
st 

constraints. The kinetic energy, T, may be expressed 
in terms of the nonholonoinic velocities and a mass 
matrix (MJ, or 

T = |{u} T IMJ{u} . (7) 

With equations (4 through 7), equations{2) and 
(3) can now be written in matrix form as 

[0] T [M]{u} + (10] T - lP])lM]{u} 


Y. = r ( v i. Y„. . Y n ,t) (13) 

where 

i - 1, 2, . . . n . 

The elements of the vector {Y} completely 
define the state of the system at any time. The 
relationships between {Y}, {y}, and t are the 
equations of motion, constraint, and auxiliary 
equations that describe the beha ior of the system. 
The elements of {Y} are initially yiven, by 
evaluating {Y} numerically, inputs into an integra- 
tion algorithm that serves to increment { Y} are 
obtained. 


= {Q'}+ lbJ T {\} (8) 

and 

lb]{q} + {« t } = {0} . (9) 

The matrix [0] is a square matrix that is, in general, 
non-singular; then, 

{<0 = UJJ-Hu} - (10) 

Now, if equation (8) is pre-multiplied by ^[/j]^ 
and equation (10) is substituted into equation (9), 
the following is obtained 

(M]{u) + ( [0] T ) ‘ (l/Sj T - [Pj) [Ml{u} 

= ((/5J T ) ' ({Q'} + lb) T {x}) (11) 

and 

(bHjl]- 1 {u} + 4 t } = {0} . (12) 


Let us define the elements of the state vector 

W 

(r T > 

{Y} = {X C } 
lT c } 

{p> 

{* T } 

U T ) 

iti 

U c > 

{ 6 Cs>J ’ (14) 

where 


Equation (11) represents n differential equations 
of motion and equation (12) represents the m 
equations of constraint. The derivation of the explicit 
forms of the matrices [P], 10], and (b) are too 
lengthy to include here and can be found in . 

Reference 2. 

STATE VECTOR 


The equations of motion and auxiliary -equations 
that characterize the docking maneuver comprise 
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and 


(y t ) 


y iit 

Y 12T 

Y 13T 

Y 21T 

Y 22T 

Y 23T 


(17) 


mass from the origin of the inertiallrame and referred 
to the inertial frame has components X^, Y and 

Z^.. The orientation of the target vehicle's body axes 


witn respect to the inertial system, is given by the 
six direction cosines, >■. . The vectors { u } , {x }, 

1J 1 L L- 


and {y } have corresponding definitions for the chase 

vehicle. The vector {p} contains six parameters that 
describe the state of the probe assembly. 


The parameters u^., v^., and are compnents 

of translational velocity of the target vehicle's center 
of mass, referred to the target body- fixed system 
(Fig. 4). The components of angular velocity erf the 
target body -fixed system are , and . 

The vector positioning the target vehicle's center of 



a 

y 



Yl 

Y2 

Yj 


( 18 ) 
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The coordinate describes the extension of 
x 

the inner cylinder (Fig. 5). The angular deforma- 
tion of the probe body resulting from elastic 
deflections of the probe support struts is described 
by the coordinates y and y . The coordinates 

y z 

(y 1 = 1, 2, 3) are the angles included between the 
pitch arms and the probe central body. 


A parameter <p. is associated with each of the 

constraint conditions. The constraints are identified 
in Table 1. 


TABLE 1. DESCRIPTION OF THE 
CONSTRAINTS (1] 



Figure 5. Probe mechanism geometry. 

The elastic deflections of the target and chase 
vehicles are expressed as a superposition of normal 
free-free elastic modes of the respective vehicles. 
The vectors } contain the normal coordinates 

1 v 

of the elastic modes of the target and chase vehicles, 
respectively. The vector {4 } contains the attitude 

L S 

control system parameters that are required to close 
the loop. 

CONSTRAINT CRITERIA 


There is a variety of combinations of constraint 
conditions that are to be implemented during the 
course of a motion simulation. The constraint criteria 
permit the mathematical model to simulate t'le 
various events that occur during the maneuvering in 
a representative manner. The motion during the 
maneuver is characterized by the bumping, sliding, 
breaking contact, and capture of the docking 
mechanisms. 


Constraint 

Parameter 

Description 

<*».. 1=1,2, 3 

jth Tension Link Constraint 
<P . = ipelo - Ipe L ; <#>j ^ 0 . 

*4 

Probe Head-on tact Constraint 
<f> 4 = -DS ■ - C ; ^>0. 

^j+4 
j = l,2,3 

jth Pitch Arm Contact Constraint 

<f>. „ = QR • m. - A ;<p. ,S0 
1+4 J J b’V4 

<P* 

Probe Head Capture Constraint 
<p, = IDS) - C/sin a ; ip, a 0 


Because of the nature of the docking system, 
there are certain constraints that cannot physically 
occur at the same time. For example, the jth 
tension link cannot be in tension when the jth pitch 
arm is in contact with the drogue lip; i. e., 4>. and 


1+4 


cannot simultaneously be equal to zero. A 


similar conflict exists between <p t and <p,. This 
means that during the motion, if <p^ = 0, then the 

corresponding constraint equation is a candidate to 
be used. When a conflict in constraint conditions 
occurs, the dominant condition is chosen. 


There is a constraint force, X, corresponding 
to each erf the constraint equations. If a particular 
constraint equation is not in effect, its corresponding 
constraint force is zero. When particular constraint 
equations are in effect, the certain elements of the 
{Y} vector are solved to satisfy the constraint 
equations. The subvector {p} is that portion of 
{*} that is committed to satisfy the constraint 
equations. 


NUMERICAL EXAMPLE 


A docking response analysis was accomplished 
considering the CSM docking to the axial port of the 
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cluster. The elastic modes of the cluster were in Figures 6 and 7. The docking force is the reso- 

included in one case and omitted in another. The lution of the constraint forces to point D (Fig. 4). 

CSM was considered rigid (inelastic) in both cases. It is noted that the maximum force resulting from 

For a given set of initial conditions, the docking the elastic simulation is significantly less than that 

force time histories were calculated and are plotted from the inelastic siir jlation. 


INITIAL CONDITION*: X - IX*. V - S IS. Z • 0 00 
SI • 110. TH - 0.00. Ft - *0 
U-1Z V • 0 w-s 



0 0.1 0.2 0.3 0.4 0.5 OA 0.7 ox a» 1.0 


FZD VS TIME 

Figure 6. Docking force time history. 


33 




W. 8. HOLLAND 


1.2 

1.0 

0.8 

0.6 

0.4 

0.2 

0 

- 0.2 

-0.4 

-<U 

- 0.8 


■1 

« 


■ 

9 

9 


9 


■1 

II 






m 


■V 





9 




Hi 

IK 




9 


i 


■1 

n 








HI 

n 






9 


1 

»■ 






9 


■ 

u 






m 


■ 

19 

■ 





9 


■i 

19 






■ 



■s 





1 


SB 


13.8 

2.0 

11.5 

8.2 ,§ 
69 iv * 2 

4, 5 * 

5 f 0-8 

2.3 3 S 

o f j" 

•4.6 

-0.4 

-8.9 

93 









m 


■ 

1 

1 

■ 

9 

9 

■ 



■ 

■ 

■ 

1 

1 

B 

i 

i 

■ 

■ 

■ 

■ 

■ 

■ 

1 

1 

■ 


■ 

■ 

■ 

■ 

■ 

9 

■ 

1 

1 


B 

■ 


9 

■ 

■ 

i 

■ 

1 

I 

■ 

i 

■ 

■ 

■ 

B 

■ 

■ 

■ 

1 

■ 

■ 

9 

9 

m 

i 

■ 

■ 

■ 

1 

n 

■ 

■ 

i 

■ 

■ 

■ 

■ 

9 

■81 

u 


'99 

■1 

■ 

99 


99 

■1 


23. » 
18.4 

13 8 3 

} 

r 

93 i 

, 

> 

o ; 

-4.6 


0.1 0-2 0.3 0.4 0.5 0.6 0.7 0.8 a9 1.0 

MXD VS TIM€ 


0.1 03 0J 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

MYOVSTIME 


ELASTIC 

INELASTIC 



40.3 

34 & 

28.8 

23.1 | 
5 

17 J % 
3 
S 

5J 1 
o 
V 

0 “ 

-IIS 

-17.3 


0 0.1 02 03 IU m Oi 17 ftl 0.1 1.0 

MZOVSTfME 

Figure 7. Docking force time history. 
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PAYLOAD ENCLOSURE NOISE REDUCTION 

By 

•3. E. Erwin 


INTRODUCTION 


During launch, aerospace structures are sub- 
jected to very high acoustic pressures on the 
exterior surface of the vehicle. These pressures 
cause the relatively thin Skinned panels and shell 
components of the structure to vibrate. Vibration 
of the exterior skin generates inter mJ noise fields 
that in turn act as fluctuating pressure sources for 
internal equipment. Internal systems such as pay- 
oads within the upper stages of the launch vehicle 
are extremely sensitive to these fluctuating pressures. 
These payloads are enclosed by cylindrical and 
conical shells, which are called payload enclosures. 
To minimize payload vibrations, it is necessary to 
maximize the noise reduction ti rough the enclosure 
structure. 

DEFINITION OF PROBLEM 


Noise reduction became a severe problem during 
design of the payload enclosure for the Skylab. When 
the Skylab program originated, the honeycomb Space- 
craft Lem Adapter (SLA) was selected as the enclosure 
structure since it was developed hardware being 
utilized on the Apollo vehicle configurations. However, 
there arose several structural problems associated 
with the SLA that involved payload interference, 
separation constraints, and a payload weight margin. 

To alleviate these problems, a skin-stringer enclosure 
was designed by MSFC. A sketch of these enclosures 
is shown in Figure 1. The SLA enclosure was 4. 3-cm 
(1.7 in.) thick honeycomb, 11.89 m (39 ft) in height, 
and 6. 61 m (21. 7 ft) in diameter. The skin-stringer 
enclosure was 12. 5 m (41 ft) in height and also 
6.61 m in diameter. The skin thickness was 0. 10 cm 
(0.04 in) and the stringer spacing was 17 p "m 
( 7 in. ) . A comparison of the noise rer n 
capabilities of these enclosures is preset.. -a in 
Figure 2. 

The SLA provided 5 to 10 dB more noise reduction 
across the frequency spectrum. It should be pointed 
out that most of the payload components and subsystems 
had been bought and/or tested to vibroacoustic levels 


uased on the SLA noise reduction. Therefore, if the 
skin-stringe” enclosure was to be used, its noise 
reduction capability must be increased to yield 
approximately the same internal acoustic levels as 
the SLA. 



SPACECRAFT LEM ADAPTER 
ENCLOSURE 



SKIN STRINGER ENCLOSURE 


Figure 1. Skylab enclosure designs (Saturn IB 
wet workshop) . 

PROPOSED SOLUTION 


Several techniques were considered for increasing 
the noise reduction of the skin-stringer enclosure. 

They are as follows: 

1. Increase enclosure mass and stiffness to reduce 
enclosure vibration levels. 
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Figure 2. Comparison of noise reduction for SLA and skin-stringer enclosure. 


2. Add absorption material to the internal sur- 
face of the enclosure to reduce the internal sound 
field. 

3. Employ a light inert gaseous medium, such 
as helium, internal to the enclosure to provide an 
impedance mismatch between the air (external) and 
the helium (internal) . 

The first of these techniques was eliminated because 
of the payload weight restriction. The addition of 
an absorptive material and the utilization of the 
helium gas were selected for further investigation. 

Helium was selected since it is the lightest, non- 
flammable gas available and therefore provides 
maximum noise reduction through a maximum imped- 
ance mismatch. The impedance of an unbounded 
acoustic medium is a function of pC (the mass density 
times the speed of sound in the medium) . The term 
pC can be thought of in terms of momentum density. 
Since C is the speed of sound in the medium or the 
speed at which mechanical energy propagates in the 
medium, pC may be regarded as the momentum per 
unit volume that may propagate through the medium. 

If a system is composed of two media having different 
values of pC, an impedance mismatch occurs at the 


interface and the inefficiency of the system is 
increased. That is, the medium with the smaller 
value of pC is capable of absorbing and transmitting 
a smaller amount of momentum than the medium with 
a greater value of pC. The excess momentum remains 
in the medium of greater pC in the form of a reflected 
wave {1}. The ratio of the impedance (pC) of helium 
to air is approximately 0. 407, and since the sound 
pressure level varies proportionally with the imped- 
ance, the expected noise reduction was predicted as 
follows: 

N.R. = 20 log - gjf - f. 6 = 8dB 
pC Air 

A survey indicated that accomplishing noise reduction 
in this manner was without precedent and would 
require careful laboratory testing. 

TEST PROCEDURE 


An applied research test program was formulated 
to evaluate and verify the noise reduction capability 
of the helium and absorptive material. A series of 
acoustic tests was conducted in the MSFC acoustic 
research facility utilizing a monocoque stainless steel 
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cylinder 91, 44 era (36 in. ( in diameter, 91, 44 cm 
(36 in. ) long, and 0, 046 era (0,013 in. ) thick and 
capped at both ends. A typical test setup is presented 
in Figure 3. The test specimen was suspended from 
the ceiling of the reverberation room by bungee cords 
and was subjected to broadband random acoustic 
excitation (2], 

Five microphones were distributed external to 
the cylinder, and five microphones were distributed 
internal to the cyiimlar. The helium purge inlet and 


outlet flow attachments were installed in the bottom 
cap. Samples taken from the bottom of the cylinder 
indicated a helium purity of at least 93 percent for 
each test run. The tests were conducted for the 
following configurations: 

1 . Unmodified cylinder with ambient air inside 
the cylinder, 

2, Unmodified cylinder with the internal volume 
purged with helium. 



Figure 3. Typical test setup. 
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3. Cylinder modified with 2. 54-cm (1-in. ) 
thick polyurethane foam applied to all interior 
surfaces and ambient air inside the cylinder. 

4, Cylinder modified with 2. 54-cm (1-in. ) 
thick polyurethane foam applied to all interior sur- 
faces and the internal volume purged with helium. 

DATA REDUCTION 


Data from the microphones were recorded on 
magnetic tape. One-third octave band analyses were 
subsequently performed on the recorded data. The 
data were plotted in decibels at one-third octave 
band center frequencies. To obtain an estimate of 
noise reduction for the cylinder, the average internal 
microphone levels were subtracted from the average 
external levels, and the resulting noise reduction 


curves were plotted. This is the form in which the 
results are presented in this paper. 

TEST RESULTS 


The results from these tests show that the use 
of helium was quite effective in providing additional 
noise reduction. These effects are shown in Figure 4 
by comparing the noise reduction measured when the 
cylinder wt.s filled with helium and when filled with 
air. The air volume critical frequency at 182 Hz and 
the helium volume critical frequency at 530 Hz are 
volume resonance phenomena that are a function of 
the media. These resonances, as expected, are 
indicated by a decrease in the measured noise 
reduction. A significant decrease in noise reduction 
also occurred near the ring frequency, which was 
calculated to be approximately 1800 Hz. It is a 



Figure 4. Noise reduction of a monocoque cylinder in a reverberant acoustic field. 
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well-known property of cylindrical structures 
that response is highest at the ring frequency, thus 
causing another decrease in noise reduction for 
both media. A large noise reduction was obtained 
in the frequency range between the air and helium 
volume critical frequencies. Above the helium 
volume critical frequency, approximately 8 to 10 dB 
additional noise reduction was obtained. In scaling 
from the test cylinder to the 6. 61 -m (21. 7-fl) 
diameter payload enclosure, the air and helium 
volume critical frequencies would shift down in 
frequency to 25 Hz and 70 Hz, respectively. Thus, 
it can be assumed that the helium will provide 
8 to 10 dB added noise reduction for the payload 
enclosure in the frequency spectrum above 25 Hz. 

The noise reduction furnished by the absorptive 
material was frequency-dependent as shown in 
Figure 5. In the frequency band from 40 to 100 Hz, 
the absorptive material provided practically no noise 
reduction above that afforded by the cylinder. 


However, from 600 to 3000 Hz, reductions of 2 to 6dB 
in sound pressure levels were obtained. An unexpect- 
ed result was the magnitude of added noise reduction 
indicated from 100 to 600 Hz. This effect was 
attributed to the added mass and stiffness caused by 
the presence of the foam since the cylinder was 
monocoque and only 0. 046 cm (0. 018 in . ) thick. 

Another series of tests was run with the 
cylinder modified with the 2.54-cm (1-in.) layer 
of foam and purged with helium. These results, 
presented in Figure 6, show that the foam and 
helium furnished approximately 25 aB in the fre- 
quency range above the helium volume critical 
frequency. These results indicate a greater noise 
reduction than the added results of the helium and 
absorptive test taken separately. This ob tion 
implies that the absorption material was mo,.- 
efficient in the helium environment. Howevtr, 
no theoretical basis was found to support tiiis 
observation. 
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Figure 6. Noise reduction of a monocoque cylinder in a reverberant acoustic field. 


CONCLUSIONS 

An added noise reduction of 8 to 10 dB can be 
obtained for a cylindrical structure with an internal 
helium environment. The noise reduction provided 
by polyurethane foam absorptive material is frequency- 
dependent. Seductions of 2 to o dB were obtained 
in the f?equency region of 600 to 3000 Hz. The 
conclusions from the helium test were the basis 
for an operational requirement to purge the payload 
enclosure with helium prior to launch. This require- 
ment specified that a 93-percent helium environment 
be maintained up to 15 sec after launch. 


Since the research effort reported in this 
paper, several additional acoustic tests have been 
oerf .- med to investigate the effect of a helium 
puige on component vibration response. A portion 
of these tests were rponsored under research 
contract NAS8-2UZ60 with Wyle Laboratories. 

Another se.ies of tests involved a large 6. 71 m 
(22-ft) diameter cylinder mounted on the Mobile 
Acoustic Research Lab (MARL) test fixture located 
at the Mississippi Test Facility. The MARL was 
Iocs- ad near the test stands and exposed to four, 
stage-8ta tic-test exposures. During these exposures, 
the large cylinder was purged with helium. The results 
from tb> »e ttwf® be icported in the near future. 
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S-ll STAGE VIIRATIQN TESTING AND ANALYSIS 
FOR THE POGO PHENOMENA 


By 

H. M. 

SUMMARY 

The Dynamics Analysis Branch of the 
Astronautics Laboratory of Marshall Space Flight 
Center conducted an investigation of the longitudinal 
o- dilation problem encountered during second stage 
'urn of the Saturn V vehicle. This report describes 
the structural test and math modeling effort that was 
undertaken. It includes details of the extensive 
experimental research program performed on the 
S-Il stage vehirle in conjunction with a definitive 
math model development utilizing the test information. 

The tests were conducted using a shortened S-ll 
stage to determine dynamic characteristics of the 
thrust structure/ center engine crossbeam/ liquid 
hydroelastic system. These data were used as 
direct input in developing the math model. Compari- 
sons between results obtained from the empirical 
model and tests are mace; these show excellent 
correlation in frequency, response, and mass 
characteristics. Also, comparisons between flight 
data and the calculated overall vehicle modal 
parameters are made to further verify the adequacy 
of the structural model. 

INTRODUCTION 

Development of large space vehicles involves 
the assemblage of vast amounts of hardware. 

Structural convergence of these components , in most 
cases, produces an operational system with many 
dynamically complex characteristics. The inflight 
dynamic response of these vehicles is dependent on 
parameters related to structural design, fluid 
mechanics, propulsion mechanisms, and control 
systems. To properly analyze dynamic problems 
that occur during flight, all of these effects must 
be taken into account. 

On December 21, 1968, the first manned 
"aturn V vehicle (AS- 503) was launched. Second 
stage longitudinal oscillation problems, known as 
POGO were first recorded during this flight. The 


Lee 

severity of this problem did not become apparent, 
however, until the flight of AS-504 (Apollo 9 mission) 
on March 3, 1970. Figure 1 shows recorded 
accelerometer response in G's versus flight time 
for three areas in the aft end of the S-H stage. 

These data were taken from oscillograms of the 
AS-504 lox sump, center engine, and outboard 
engine during S-ll stage burn. 

An analysis of this flight data and some previous 
test data revealed that the center engine crossbeam 
responded at a frequency near 18 Hz throughout 
flight. The response of the beam was erf low ampli- 
tude until near the end erf second stage burn. At 
this time the system became unstable, and the large 
amplitude oscillations shown in Figure 1 were 
encountered. It was hypothesized that one of the 
predominant lox bulkhead modes increased in 
frequency with decreasing liquid level and coupled 
with the crossbeam resonance to cause thir 
instability. This would apparently indicate that as 
the lox bulkhead frequency approached the 18 Hz 
region, the total system would enter into a closed 
loop, regenerative mode, producing high response 
levels in the aft section of the S-D stage. 

Because of the impact and complexity of the 
problem, an extensive experimental program was 
initiated. Figure 2 exhibits a flow chart outline 
of the total S-H oscillation program. This paper 
presents, however, only a summary of the testing, 
modeling techniques, and data comparisons utilized 
in updating, assessing, and verifying the structural 
model. 

DESCRIPTION OF TESTS 


Dynamic testing of die S-H stage was composed 
of three phases. Figure 3 shows the structural 
differences between the three phases as well as the 
method of dynamic excitation. The approximate 
location of the accelerometer instrumentation is 
also shown for ea a case. It should be noted that for 
phase 1 testing, the aft skirt was fixed to ground, 
but during phase II and 01 testing the structure was 
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Figure 1. AS- 504 vibration amplitudes. 


supported on a low frequency air suspension system. 
The lox tank was filled with de-ionized water and 
drained in several incremental levels. Modal surveys, 
which include sinusoidal sweeps from 5 to 50 Hz and 
resonant frequency dwell tests, were then performed 
at each of the levels. 

The primary objective of the phase I sequence 
was to empirically investigate the axisymmetric 
modes of the aft lox tank bulkhead as a function of 
propellant level, corresponding to levels prior to and 
through the period of longitudinal oscillation observed 
during the AS-504 flight. The investigation included 
the acquisition of such dynamic parameters as mode 
shapes, damping frequency, generalized mass, and 
tank bottom pressures. The phase H and phase in 
test efforts were designed to determine bulkhead/ 
thrust structure/crossbeam modal interaction. 

Phase IQ also included the installation of the five lox 


feedlines capped and filled with water. Tim goal 
of these two phases was to determine the following 
properties: (1) predominant system mode shapes, 

(2) system damping, (3) modal frequencies, 

(4) generalized mass, (5) effects of the lox feedlines, 
and (6) tank bottom/lax line pressures. 

The phase 1 test proved to be an excellent setup 
for determining the axisymmetric nydroelastic 
modes of the lox tank bulkhead. Accelerometers 
were placed along several meridians over a 180- 
degree section of the bulkhead. Response data 
obtained at the same level on each meridian were then 
averaged to yield mode shapes. Figure 4 shows the 
first three bulkhead resonant frequencies versus flight 
time as related to the AS-504 vehicle S~n stage burn. 

Phasi 11 and 111 data were obtained from the 
structural intera ction of the lox tank bulkhead, aft 
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Figure 2. S-II oscillation overall program flow. 


skirt, thrust structure, center engine crossbeam, 
and several other components. Presented in 
Figure 5 are the first five longitudinal S-Q system 
modal frequencies versus AS- 504 second stage flight 
time. Interesting modal frequency patterns can be 
noted in the figure because erf the complex interaction 
of the structural components. These frequency 
versus time plots represent only one of the many 
types of forced response information obtained from 
the S-II testing program. It was from the data that 
the math model of the aft end of the S-II vehicle 
was defined. 

MATH MODELING TECHNIQUE 


A multi-degree-of -freedom model of the phase I 
test configuration was developed to simulate the 
resDonse of the structure/fluid system. The modeling 
technique used is given in Reference 1. The structure 
and fluid were assumed to be a damped structural 
system subjected to harmonic forced motion. Test 
response data are used as direct input to the equation 
of motion, from which the uncoupled mo'*e shapes 
are obtained. The corresponding mass and stiffness 
matrices can be calculated using the condition of 
orthogonality of modes. 

The phase II modeling utilized tie same technique 
as described above, incorporating not only the lax 


tank but also the thrust structure, outboard engine, 
and center engine crossbeam. The governing 
differential equation in generalized coordinates is 

{%) + + 

l<p]^ (F sin ait} , (1) 

where 

q - nth model coordinate 
n 

u> q = natural frequency of nth mode 

{ d>} = nth modal column 

{ F sin cot} = applied force at frequency u> 

S = damping of nth mode 
n 

Mg = generalized mass of nth mode, 
n 

The steady-state solution of equation (1) is 
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DAMPING MODE SHAPES 

Figure o. Test definition. 

in which 

d = tan -1 

T n 


and 

^ *Sn 
where 

F = applied force at < oint r in the nth mode 
rn 

p - modal displacement at point r in the nth mo^e 
rn 

= modal displacement at point Sin the nth mode 
Fg n = real displacement at point Sin the nth mode. 




MfUBMKISCI 


Figure 4. Phase I frequency versus time. 
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Figure 5. Phase n and in frequency 
versus time. 


This substitution then results in the following 
equation: 



Sn 

If — — is treated as forced response, one can 
r 

m 

assume that it is comprised of the summation of the 
response contributions from each mode being ana- 
lyzed. With this in mind, the following system of 
equations is obtained for response with an applied 
force at mode 1 with a forcing frequency of <y = u> 1 : 


Cn<Pu + C + . . . + C 1q ^ = ^ 
c n0u$ii + CuQiz't’ii + • • • + c ln ^ ln ^2n = Fj 


C «*«*nl + ClJ0, ^n 2 + • ‘ ‘ + C m*ln*nn “ 

The constants C represent the coefficients of 

the modal displacement products <4 <p seen in 

rn sn 


equation (3). The subscripts i and j refer to the 
forcing frequency w. and the modal resonant frequency 


ur, respectively. 


C. 


can be more specifically 


defined as follows: 


C.. = 

ij 


sin (w.t +>!■.) 

rr* — L 




(4) 


A similar set of equations may be obtained for every 

forcing frequency <u - u>. to u> = w . 

l n 


Careful examination of equation (4) reveals 
that the modal displacements are the only unknowns. 
All other parameters can be obtained for the test 

data. These unknowns [" d» <t>„ 1 supply the informa- 
L rm Sn J 

tion needed to define independent normal modes [<f>J 
for the system. 


Utilizing the conditions of orthogonality, the 
mass and stiffness matrices can be found in the 
following manner: 

[M] =[« T ] _ n ‘ [XJMn' (5 

IK1 = [*T D Mg nJHn ' (6 

Several fill conditions were selected from each 
test phase. Modal frequencies, normal displace- 
ments, damping, and generalized mass from the 
dwell test data were input into equation (1) , and 
normal modes were then substituted into equations 
(5) and (6). 


COMPARISON OF TEST AND MODEL DATA 


One check made on the developed model was a 
comparison of calculated response versus frequency 
with frequency response data obtained from the test. 
The model was derived using data obtained at reso- 
nant dwell points. As expected, exact agreement 
was shown at those points, as seen in Figu- 6 and 7 
for phase I and phase n comparisons, respectively. 
Also good agreement is shown at points other Uian 
resonance through the entire frequency range of 
interest. This indicates that the model adequately 
simulated the hydroelastic test system. 

The total longitudinal mass of the system was 
determined from the generated mass matrices and 
is shown in Table 1 for the Phase II conditions. The 
total mass check was important, since this model 
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Figure 6. Phase I — model/ test comparison 
(205-sec flight time) . 


would eventually be incorporated into a flight 
model of the Saturn V vehicle. No restraint had 
been placed on the magnitude of the masr matrix, 
but the table shows good agreement with actual 
test vehicle weights. 

COMPARISON OF VEHICLE ANALYSIS 
AND FLIGHT DATA 

It should be noted that the phase III test data were 
not used in deriving a final model for use in vehicle 
analysis. A comparison between the phase II and 
phase III data did not reveal an anomaly induced by 
the addition of engine fee.Uines. It was, therefore, 
concluded that the line definition that had been 
previously used would be ('.corpora ted and that 
phase II data adequately defined the thrust structure/ 
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Figure 7. Phase II — model/test comparison (205-sec flight time). 


TABLE 1. COMPARISON OF ACTUAL WEIGHTS WITH WEIGHTS CALCULATED FROM MASS MATRICES 


S-II Flight 
Time, sec 

Actual 

Weight, kg (lb) 

Mass Matrix 
Weight, kg (lb) 

Percent Error 

205 

207 445 (457 330) 

218 051 (480 711) 

+5.11 

313 

97 447 (214 830) 

119 730 (264 000) 

+22.8 

328 

84 001 (185 187) 

82 826 (182 597) 

-1.39 

341 

72 489 (159 588) 

71 851 (155 977) 

-2.26 

352 

62 882 (138 629) 

62 847 (138 332) 

-0.02 

361 

55 772 (122 954) 

55 903 (123 243) 

+0.02 

361 

52 853 (116 496) 

52 759 (116 311) 

-0.01 

1 370 

50 197 (110 663) 

51 482 (113 496) 

+2.56 
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center engine crossbeam/ bulkhead system. A model 
derived from these data was then incorporated into the 
total Saturn V vehicle (Fig. 8), which contains spring- 
mass representations of the remaining flight structure. 
Figure 9 shows a comparison of a frequency analysis 
of the Saturn V S-II stage burn model versus AS-504, 
AS-507, and AS-508 flight data. The first three 
modes seen on the chart are major vehicle longitudinal 
modes. The next three modes represent the first and 
second predominant bulkhead modes and the predomi- 
nant crossbeam mode, respectively. It is apparent 
that the flight data generally follow the predominant 
second bulkhead mode. However, no definite conclu- 
sion can be drawn from this analysis , since ■•h 
parameters as forcing function, control systems, 
and engine transfer functions were not incorporated. 

Approximate flight time for the POGO instability 
on AS-504 was at 340 seconds into the S-II stage burn. 



The AS-507 flight did not experience this instability 
because of the center engine shutdown at about 
60 seconds, before the end of S-II burn. This was 
done to open the unstable rege nerative loop 
developed on flights such as AS-504. The AS-508 
flight also had a scheduled early center engine 
shutdown; however, POGO dee-eloped at about 170 
seconds into S-II burn, and an unscheduled center 
engine shutdown was effected 1 ecause of excessive 
vibration leveis. 


CONCLUSIONS 


It has been shown that the aft end of the S-II 
stage, including such complex components as 
bulkhead/ fluid, thrust structure engines, anu cross- 
beam, was successfully tested and mathematically 
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Figure 8. S-II flight m< del. 
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F ; g'ure 9. Flight data versus flight model. 


modeled. The results yielded an excellent frequency 
response and mass characteristics correlation. 

These models were incorporated into overall flight 
vehicle models and analyzed for modal parameters. 
Preliminary comparisons of frequency analysis with 
actual flight data indicate that the total structural 
modeling program has produced a good representation 
of the flight vehicle, -dual POGO instability 
analyses will be performed later utilizing these test- 
developed models along with other propulsion and 
control system transfer function modeling data. 

Future Saturn V flight vehicles will also be evaluated 
for this oscillatory problem when the model is able 
to reproduce all responses observed on previous 
vehicles. 
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